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We investigated binding of hydrogen atoms to small Polycyclic Aromatic Hydrocarbons (PAHs) - i. e. graphene 
dots with hydrogen-terminated edges - using density functional theory and correlated wavefunction techniques. 
We considered a number of PAHs with 3 to 7 hexagonal rings and computed binding energies for most of 
the symmetry unique sites, along with the minimum energy paths for significant cases. The chosen PAHs are 
small enough to not present radical character at their edges, yet show a clear preference for adsorption at the 
edge sites which can be attributed to electronic effects. We show how the results, as obtained at different level 
of theory, can be rationalized in detail with the help of few simple concepts derivable from a tight-binding 
model of the tt electrons. 
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I. INTRODUCTION 

Graphene, the recently discovered two-dimensional 
form of carbon^, is a promising material for a future 
carbon-based nanoelectronics. Its peculiar tt — tt* elec- 
tronic band structure, with a linear energy dispersion 
close to the Fermi level, introduces subtle quantum 
pseudo-relativistic effects in the low-energy charge car- 
rier dynamics which hugely impact on the transport 
propertiea^'Hl. This results, e.g.., in a robust anoma- 
lous quantum Hall effect^ ^, a universal conductivity 
minimum^ and ballistic transport which can reach the 
micrometer scale^. From a practical point of view, the 
substrate thickness, the high mobility of its charge car- 
riers and their (high-field) high saturation velocity rep- 
resent attractive features for the chip-makers. Nanos- 
tructuring, however, is needed for applications, e.g. for 
devising graphene-based logic transistors where a band- 
gap is needed to achieve high operational on-off ratios. 
Graphene Nanoribbons (GNRs) can be cut which show 
either semiconducting or metallic properties, the latter 
coming with edge states of unusual magnetic proper- 
ties, possibly leading to carbon based nanomagnets^ESl. 
Likewise, Graphene Dots (GDs) can be designed to 
have specific electronic structures and transport prop- 
erties, by acting just on their shape and their con- 
nectivity. GDs hav e bee n suggested for realizing spin 
qubits-^, spin filters^^^ and spin- logic devices^^, and 
proposed as biomedical imaging agents^^ and light ab- 
sorbers for photovoltaic^^. Transport properties have 
been measured on a variety of dot devices carved en- 
tirely from graphene by high-resolution electron-beam 
lithograph}^^''. 

Most of these properties arise entirely from the tt elec- 
trons and remain unaltered when saturation of the dan- 
gling <j bonds occurs, e.g. in forming Polycyclic Aromatic 
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Hydrocarbons (PAH). The latter offer an enhanced chem- 
ical stability, and their nanostructuring (energy level ar- 
rangement, interfacing with other materials, etc.) can be 
realized with the help of well- developed carbon chemistry 
methods. They have been used as building blocks for 
atomically-precise nanoribbon fabrication^^ and, in prin- 
ciple, may form the basis for a bottom- up approach to 
realize arbitrarily complex carbon-nanostructures. PAHs 
have also been investigated in many other fields, from 
petroleum chemistry to astrochemistry. For instance, in 
the interstellar medium (ISM), i.e. the extremely rar- 
efied medium which fills the space between stars, the 
observed abundance of molecular hydrogen cannot be 
explained by direct gas-phase routes involving H atoms 
only, rather is believed to occur on the carbonaceous sur- 
face of dust grain d-^^'^^' and small carbonaceous particles. 
PAHs, which are estimated to lock up ca. 15% of the 
interstellar carbon, have been suggested as possible cat- 
alysts for H2 format ioiP^MJ^ 

In this work, we investigate the reaction of atomic hy- 
drogen with a nu mber of PAHs, complementing previ- 
ous related studies -* ^^ * ^^" ^ which showed preference for 
addition at the edges of selected PAH molecules. The 
main aim of this study was to emphasize the importance 
of substrate relaxation ("geometrical") effects in deter- 
mining a preference towards the edges. To this end, we 
selected substrate PAH molecules with relatively small 
(sub-nanometer) dimensions, in such a way to prevent 
any enhanced chemical reactivity at the edges due to a 
true radical character (single occupation of a semilocal- 
ized edge state), as it occurs for instance at the edges 
of wide zig-zag GNRs. However, as we shall see in the 
following, some edge localization is always present. This 
provides an enhancement of the edge reactivity which is 
of purely electronic origin and can be easily understood 
in terms of few concepts derivable from a tight-binding 
(Hiickel) model for the tt electrons. 

In addition, depending on the number of carbon atoms 
available for the tt electron system and their connectivity, 
the systems considered can also show a marked suhldXtice 




Figure 1. Left: Two- (red) and three- (blue) coordinated edge 
sites in the benzo[ghi]perylene molecule {E and F sites in the 
main text). Right: E sites having different hypercoordina- 
tion number (as indicated), along with their hypercoordinated 
partners (white circles). 



preference due to the "alternating paths" followed by (un- 
paired) itinerant electrons in graphenes {i.e. due to the 
presence of st aggered midgap states). This is similar to 
graphen d^^ l ^^ l, where these states form the basis for a 
preferential sticking mechanism^^ ^^ forming para-dimers 
{i.e. two H atoms on opposite corners of the same ring). 
We therefore distinguish two classes of PAHs according 
to whether the number of sites in each sublattice is bal- 
anced or not, and show how a final set of rules governing 
the site reactivity results from the interplay of different 
electronic effects. The basic concepts underlying these 
rules equally apply to larger systems, and thus allow one 
to easily predict the chemical reactivity of sp'^ carbon 
nanostructures with monovalent species forming covalent 
bonds with the substrate. 

Importantly, in the present study we also take advan- 
tage of the modest size of the systems investigated, and 
exploit the unique opportunity of assessing the quality 
of the results of commonly used Density Functional The- 
ory (DFT) methods in investigating chemically-derived 
graphene structures. This is done here by complement- 
ing the DFT data with those obtained by using more ac- 
curate correlated wavefunction techniques. En passant, 
we briefly discuss the magnetic properties of pristine and 
hydrogenated PAHs, which turn out to be well predicted 
by Lieb's theorem^^, in agreement with previous studies 
on triangularly and hexagonal ly sh aped GDs^- and other 
defective graphenic structure u^^ l ^^ l 

The paper is organized as follows. Section I introduces 
some basic properties of 7r-conjugated carbon systems 
which underlie the presentation of the results given in 
Section III, after Section II has provided the computa- 
tional details of our calculations. Section IV summarizes 
and concludes. 

Notice that in the following we adopt a surface science 
terminology, whereby "adsorption to the substrate" (here 
meant to be chemisorption) is used interchangeably with 
"binding to the molecule". 



II. BASIC PROPERTIES OF tt ELECTRONS IN sp'^ 
CARBON STRUCTURES 

Carbon sp'^ structures like graphene, GNRs and GDs, 
are characterized by a bipartite lattice where two distinct 
sublattices, A and B, can be identified such that each A 
site is connected to B sites only and viceversa. This has 
important consequences in the one-electron spectrum if, 
as it is the case for such structures, the transfer (hopping) 
energies beyond the nearest-neighbors are of secondary 
importance and the orbital overlap can be neglected. Un- 
der such circumstances, indeed, it is not difficult to prove 
that the tight-binding Hamiltonian for the Pz orbitals 
of the TT electron system has a simple symmetry. Such 
Hamiltonian reads as 

H^^ = Yl ^^A^J + ^'^' = ^^B + Hab 

where a^(aj) annihilates (creates) an electron in site i 
of the sublattice A (similarly for bj{bj) and B sublat- 
tice sites), tij is the hopping between sites i and j, and 
the on-site energy (the energy of an isolated Pz orbital) 
has been set to zero. Bipartism is responsible for its 
(off-) block structure - as emphasized here with the in- 
troduction of Hab and Hba which collectively describe 
the transitions A ^ B and B ^ A, respectively - and 
easily leads to a symmetric spectrum around e = 0. The 
latter is also the position of the Fermi level with one elec- 
tron per site (half- filling) , and for this reason the above 
symmetry is also called electron-hole symmetry. In con- 
junction with the spatial symmetry, the presence of such 
symmetry is at the origin of the conically shaped band 



structure of graphene close to the Fermi level^, with in- 
terestmg consequences on band-engmeermg"^'^ "^^ and on 
the chemical reactivit}'^^. Here, to clarify the connection 
with chemical reactivity, we focus on some simple results 
on the shape of low energy {i.e. close to the Fermi level) 
orbitals that directly follow from such electron-hole sym- 
metry. 

Edge localization and hypercoordination. Low 
energy orbitals show a marked tendency to localize on 
edge sites, as can be easily seen at the tight-binding le vel. 
To this end, we perform a lattice "renormalization'' ^^1^ 
and focus on one sublattice only (say A) and on the 
"renormalized" Hamiltonian H = Hab Hba- The renor- 
malized energies e^ are simply related to the eigenvalues^ 
ef of H^^ , ef = ±>/ei7 and the renormalized lattice 
is a triangular lattice (the sublattice A of the original 
system) with hopping t^ [assuming tij = t for simplic- 
ity] and on-site energies t^Z^, where Zi is the coordi- 
nation number of the i — th A site in the original lat- 
tice. An edge necessarily has undercoordinated {Z = 2) 
sites, hence the ground- state of the renormalized lattice 
{i.e. the highest occupied/lowest unoccupied molecu- 
lar orbital [HOMO/LUMO] pair of the original lattice) 
naturally tends to localize on these sites which present 
the lowest on-site energy. In the following, we name 
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Figure 2. Balanced (a-d) and imbalanced (e-g) PAHs investigated in this work, shown via one of their possible Lewis 
structures (carbon atoms are at the vertexes of the hexagons, and are meant to be saturated with hydrogen atoms -not shown- 
if undercoordinated). (a) Pyrene, (b) anthanthrene, (c) benzo[ghi]perylene, (d) coronene, (e) phenalene (peri-naphtene), (f) 
benzo[c]pyrene and (g) benzo[c] anthanthrene. Also indicated a labeling systems for the adsorption sites considered in this work, 
E and G for "edge" and "graphitic" sites, respectively. A prime is used for edge sites with hypercoodination number ^ = 2 and a 
star is used in (e-g) for the majority sites, either E or G, where the unpaired electron (dot) localizes. See Section ni] for details. 



E these two-coordinated edge sites, to distinguish them 
from those three-coordinated sites which are also present 
at an edge {F sites), see Fig. fl] Importantly, we expect 
that low-energy orbitals localize on E sites and, among 
these, on those sites which show the largest number of un- 
dercoordinated neighbors in the renormalized lattice (or, 
equivalently, next-to-nearest E neighbors in the original 
lattice) to hybridize with. As is shown in the following, 
this latter number turns out to be an important parame- 
ter ruling the reactivity of the edge sites; for this reason 
we call it the hypercoodination number (^). Fig. [T] right 
panel, reports some illustrative cases. 

Midgap states and spin alignment. Obviously, 
energy levels at e = 0, if present, play a major role at 
half-filling in determining the reactivity and the mag- 
netic properties of the GDs. It is instructive to see when 
this situation occurs, as this also adds further constraints 
on the spatial behaviour of the low energy orbitals. In 
general, the number of these "midgap" states is deter- 
mined by the site-connectivity but their occupancy (spin- 
alignment) is solely determined by the sublattice imbal- 
ance. This follows from a rigorous result proved by Lieb^^ 
for the realistic (repulsive) Hubbard model having H^^ 
above as one-electron Hamiltonian: Lieb's theorem states 
that at half-filling the ground-state spin S is given by 
S = \nA — ^s|/2 where ua and ub are the number of 
sites in sublattice A and 5, respectively. 

Typically, the number of midgap states matches the 
sublattice imbalance, since this is enough to allow for 
I^A — ^b] linearly independent eigenvectors of H^^ at 
zero energy, all with null amplitudes on the minority 
sublattice sites^. [Accordingly, in this case, Lieb's theo- 
rem above becomes a sort of Hund's rule applied to the 
midgap states.]. This a simple algebraic result: for, let 



^A > ^B and lip) = Xli<^i \^i) t>e a trial solution (here 
\ai) = a] \0)). At zero energy, ^. {bj\H\ai) ai = must 
hold for j = 1, ..n^, which is a set of ub equations for the 
riA > TLB unknowns a^ having (at least) ua — tlb linearly 
independent solutions. This also shows that V^'s localize 
on the A lattice sites. 

More generally, the concept of non- adjacent sites in 
a A^-site bipartite system helps counting the number of 
midgap states^^. We say that two sites are non-adjacent 
if they are not bound (connected by a transfer integral) 
to each other; for instance, two sites on the same sub- 
lattice are non-adjacent. Clearly, there exists a maximal 
set of non-adjacent sites and we call a the sites in this 
set, and /3 the remaining ones {na-, rif^ = N — n^ ii^ num- 
ber, respectively). Each site a binds at least one site /3, 
otherwise it would represent a completely isolated site. 
Arranging one electron per site a, however, we can form 
at most 71/3 bonds at a time, and therefore we are left 
with T] = fia — rif^ = 2na — N unpaired electrons. Equiv- 
alently, we end up with r] midgap states localized on the 
maximal set of non- adjacent sites. The case of a sublat- 
tice imbalance discussed above is a special result of this 
rule which, as is evident from the discussion above, can 
be equivalently re-phrased by defining r] to be the num- 
ber of unpaired electrons in the Lewis structure (s) with 
the maximum number of tt {i.e. double) bonds. 

We thus see that, in addition to the edge localization 
discussed above, depending on the number of sites and 
their connectivity, there may exist topological constraints 
which force the carbon sp^-system to have zero energy 
states. The latter localize on specific lattice positions 
which are easily identifiable by inspection. 

The systems. In the following we mainly focus on 
structures where the sublattice imbalance is the only 
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Figure 3. Convergence tests on the active space used in the 
MCQDPT calculations. The energies for H atom adsorption 
are reported for different sites as functions of the number n 
of active electrons in a (n,n) correlation scheme. Left: grey, 
red and blue symbols for sites Ei , E2 and G of the phenalene 
molecule (structure (e) in Figpl. Right: grey, red, green and 
blue symbols for sites Ei, E2, E3 and G of pyrene (structure 
(a) in Fig|2|. Horizontal lines mark the values obtained at 
the DFT level. 



source of migdap states, and call them balanced {S = 0) 
or imbalanced {S > 0), accordingly; in particular, only 
structures with one unit of imbalance are considered, i.e. 
they all have S = 1/2, as suggested by Lieb's theorem 
and confirmed by our calculations. The considered struc- 
tures are shown in Fig|2J together with a labeling system 
for the sites investigated, which distinguishes the (two- 
coordinated) edge sites from the graphitic sites, E and G 
in Fig|2j Sites at the edges which are three-coordinated 
(F) are in between the two categories and will not be 
considered in the following. With this exception, all the 
symmetry unique sites were investigated for binding of 
a H atom, with the methods described in the following 
Section. 

Notice that Fig |2] further distinguishes those edge sites 
which have the largest possible hypercoordination num- 
ber (^ = 2) with a prime and, where appropriate, 
identifies with a star the majority sites (either edge or 
graphitic) where the midgap states are expected to lo- 
calize. As is shown in the following these labels help 
identifying the sites with the highest hydrogen affinity 
{i.e. the sites with the largest binding energy and the 
smallest barrier to binding). 



III. COMPUTATIONAL METHODS 

For each of the selected PAH molecules we computed 
the binding energy of a hydrogen atom to the sites labeled 
in Fig J2] according to 

Ebind = FpAH + % - ^PAH-H 

with two different electronic structure methods. PAH 
structures were optimized at the (unrestricted) Den- 
sity Functional Theory (DFT) level using the popu- 
lar B3LYP hybrid exchange-correlation functional with 
Dunning's double- valence, atom-centered basis set of the 



correlation-consistent type (cc-pVDZ), as implemented 
in GAUSSIAN 03^^ On the DFT-optimized structures 
single-point wavefunction calculations were performed 
with the same basis-set. These are of the multi-state, 
multi-reference perturbation theory type according to 



the scheme of HiracPHEl ^^^ NakancPEH called Multi- 
Configuration Quasi-Degenerate Perturbation Theory 
(MCQDPT) and implemented in GAMESS'^^ In this 
scheme dynamical correlation is introduced in a Multi- 
Configurational (MC) wavefunction by properly defin- 
ing a reference one-electron Hamiltonian based on this 
wavefunction and computing the second-order pertur- 
bation correction. The chosen MC reference wavefunc- 
tion was of the Complete Active Space Self- Consistent- 
Field (CASSCF) type, where n valence electrons are dis- 
tributed in m orbitals ( CAS (n,m) in the following) and 
self-consistency is reached in a variational optimization. 
In principle, for the PAHs above a consistent procedure 
would require to put all the tt electrons of the substrate 
molecules and that of the H atoms in the same num- 
ber of orbitals. This is of course impracticable for all 
but the smallest molecules, and we therefore resorted to 
an orbital localization procedure which takes advantage 
of the local character of the bond formation process. We 
started from Pipek-Mezey ROHF localized orbital^^ and 
included in the active space the a orbital describing the 
formation of the C — H bond and the tt orbitals localized 
on the sites which are nearest neighbors of the binding 
site. This gives rise to typical CAS(9,9) or CAS(8,8) 
MCSCF wavefunctions and active spaces for the pertur- 
bation correction. For the smaller PAHs, we performed 
some convergence tests on the size of the active space, 
see FigJ3] for an example. Finally, we also performed 
plane-wave based, periodic DFT calculations with the 
help of the VASP code^^^^, wit h para meters similar to 
those used in our previous workP^I^ but adapted to a 
cluster calculation. Briefly, we adopted a 20 Ax20 Ax20 
A cell and a 700 eV energy cutoff, with a 1x1x1 T centered 
/c-point grid. Inner electrons were frozen by the projec- 
tor augmented wave^^ ^^ (PAW) approach, and exchange- 
correlation effects were handled with the Perdew-Burke- 
Eznerhof^^ (PBE) functional in its spin polarized version. 



IV. RESULTS AND DISCUSSION 

A. Graphitic vs. Edge sites 

We start by showing the preference for adsorption 
on the edge s ites which was already noted by several 
author j^ ^ l ^^ l ^^" ^. Fig. [i] shows the computed binding 
energy for all the E and G sites of the structures (a-d) 
of Fig. [2J as obtained in the ground-state spin manifold 
of the total systenPl, S = 1/2. The DFT results (black 
histograms) compare very well with the available litera- 
ture data. For instance, for the pyrene molecule we find 
1.53, 1.67 and 1.09 eV for the sites ^i, E2 and Es which 
compare well with the values 1.50, 1.61 and 1.06 eV re- 
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Figure 4. Binding energies for E (left) and G (right) sites in the structures (a-d) of Fig. [2] DFT and MCQDPT results are 
represented as black and green histograms, respectively, according to the labeling system of Fig[2] 
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Figure 5. (a) Optimized structures for the singly hydro- 
genated coronene molecule. Left and right panel for adsorp- 
tion on the G and E site, respectively, (b) Reorganization 
energy for adsorption of a H atom in the indicated sites of 
the benzo[ghi]perlyene molecule. 



cently obtained by Rasmussen et alP^ with a real-space 
implementation of the DFT-PBE level of theory. 



Clearly, a striking difference between E and G sites 
is apparent from Fig. |4] binding energies at an E site 
can be as large as twice the binding energy for a G site. 
The latter, on the other hand, compare rather well with 
the valu e of th e hydrogen atom adsorption energy in 
graphene^^lSSI^ and graphit^^^^. This simple finding, 
together with a corresponding behaviour for barrier en- 
ergies to be discussed below, already suggests that the 
edges of realistic samples could be active sites where hy- 
drogenation starts and propagates into the bulk. 



B. Geometric vs. electronic effects 

Before analysing the results in details, we show here 
that "geometrical" effects per se cannot explain the dif- 
ferent behaviour of edge and inner sites evident in Figj4] 
Binding of a H atom on a sp^— carbon atom requires a 
sp'^ -^ sp^ rehybridization which leads to a tetrahedral 
reorganization of the bonding partners, as is shown in 
Fig. [51(a) for the case of the coronene molecule. Without 
such re-arrangement of the local environment no binding 
would occur: a local substrate relaxation is essential to 
"prepare" the electronic structure for binding, but this 
too is affected by the overall electronic structure which 
is always dominated by molecular orbitals spreading all 
over the molecule. 

A simple (but wrong) argument would suggest that the 
same local re- arrangement which occurs upon bonding 
(but without the H "probe") requires less energy for an 
edge than for an inner carbon atom, since in the first 
case at least one of the bonding partners is a monovalent 
species not embedded in the molecular network. We can 
define this reorganization energy as 

Er = E:^{PAH)-E,,{PAH) 

where Eeq{PAH) is the energy of the pristine molecule in 
the equilibrium configuration and E*q{PAH) is the en- 
ergy of the molecule in the same distorted configuration 
that it takes when binding the H atom. In contrast to 
the expectation above, we find that Er for an E site is 
always larger than that for a G site. For coronene, for 
instance, we obtain 1.40 eV and 1.04 eV, respectively, at 
the DFT level of theory, and similar values are found for 
all the structures considered in this work: the reorgani- 
zation energy is '^1.4±0.1 eV for E sites and ~1.0±0.1 
eV for G sites, see for instance Fig. Islb) for the case of 
the benzo[ghi]perylene. We thus see that the preference 
in binding a H atom to an edge site occurs despite the 
larger reorganization energy needed at these kind of sites. 
This allows us to conclude that this preference is due to 
the electronic effects introduced in Section HH 
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Figure 6. Binding energies as functions of the site normalized 
populations of the substrate HOMO. Black, red, green and 
blue symbols for structures (a-d), respectively. Also indicated 
the hypercoordination number of the edge sites. 



C. Hypercoordination 

Next we discuss the results of Fig. |4] in detail since, 
apart from the overall behaviour, the binding energies 
can take quite different values depending on the site they 
refer to. A closer inspection reveals that the values for 
the interesting E sites correlate very well with the hyper- 
coordination number introduced in Section [TT| the lager 
is the hypercoordination the larger is the binding energy. 
This can be made evident by reporting the results of Fig. 
[4] as functions of the site populations Pi of the HOMO; 
the latter are meant here per spin species, and were ob- 
tained by a Mulliken analysis of the molecular orbitals 
of the pristine molecules, as computed with a restricted 
Kohn-Sham determinant. This is shown in FigJ6] where 
the populations have been normalized to the values they 
would have if the HOMOs spread over all carbon atoms 
(1/A^). Fig. |6] shows that the binding energies correlate 
well with Npi. The trend is roughly linear and differ- 
ent for the E and the G sites, but we did not attempt 
to extract any behaviour because of the limited num- 
ber of data available. More importantly. Fig. [6] shows 
that the energies correlate well with the hypercoordina- 
tion number ^ of the site, particularly if the comparison 
is made between sites of the same molecule. As already 
emphasized above, this number can be readily obtained 
by simply inspecting the carbon structure under study. 



D. Imbalanced structures 

Next we move to the more complicated situation (the 
doublet structures (e-g) of FigJ2|) where topological con- 
straints lead to the appearance of zero-energy modes and 
additional "localization". Analogously to the results of 
Fig. [4J we find also in this case a clear distinction be- 
tween edge and graphitic sites. This is evident from Fig. 
[7| where we report the binding energies for the structures 



(e-g) on a larger energy scale than the one used in Fig. 
|4] This is one of the consequences of the additional elec- 
tronic effect due to the appearance of the (singly occu- 
pied) midgap state: binding of two radical species only 
requires coupling of their unpaired electrons and is thus 
typically much more energetic than in the case where a 
bond has to be broken. A further consequence is a rough 
splitting of the results into two "branches", according to 
whether the relevant site belongs or not to the majority 
set (red and blue blocks of results in Fig. [t]). Notice that, 
according to Lieb's theorem, in the first case the result- 
ing total spin state is a singlet, whereas in the second 
case is a triplet. This indeed what we find: Fig J7| shows 
that for adsorption of a H atom on a minority site the 
binding energy in the triplet state is larger than in the 
singlet. Not shown in the figure, we also checked that 
adsorption on a majority site occurs more favourably in 
the singlet manifold; this is true for all cases considered 
but the site Gi of structure (g) where we find that H 
binds more favourably in the triplet state^^. 

We thus see that, in the case considered in this section, 
the energy ordering arises from the complicated interplay 
between coordination, hypercoordination and topologi- 
cal frustration. For this reason, in plotting the results 
as functions of the normalized populations, analogously 
to Fig. [6) we consider separately the majority and the 
minority sites, reported in the left and right panels of 
Fig. [8J respectively. We see now that a good correlation 
between the binding energies and the HOMO popula- 
tions is found only for the majority sites, nevertheless 
the hypercoordination number remains a good parame- 
ter for establishing the right energy ordering within each 
category: the binding energy is found to monotonically 
increase when increasing ^. In general, majority sites 
show larger binding energies of minority sites with the 
same coordination number {i.e. either E or G) but, even 
for the same molecule, a large hypercoordination may 
offset the topological frustration of a minority site. For 
instance, the (minority) site Ei in structure (g) shows a 
larger binding energy than the (majority) site Es; notice 
though that the "expected" ordering is restored if com- 
parison is made between results for the same spin man- 
ifold. In general, however, the most favoured (relevant) 
final hydrogenated structures are always easily identified: 
they are obtained by binding a H atom to the majority 
E sites with the largest hypercoordination number. 

Notice further that imbalanced structures also arise af- 
ter a H atom has been adsorbed onto any of the balanced 
structures (a-d), since formation of a CH bond effectively 
removes one carbon pz orbital from the tt network and 
thus acts as a vacancy. In this case, hydrogen bonding 
to form a dimer follows the same rules. For instance, 
Rauls and Hornekaei^^ used DFT-PW91 to systemati- 
cally investigate hydrogenation of coronene up to satura- 
tion. They found that addition of a H atom to the most 
stable H-coronene structure {i.e. with a first H bound 
to di E site) is most favoured in the ortho-edge position, 
i.e. on the E site which is nearest neighbor to the first 
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Figure 7. Binding energies for E (left) and G (right) sites in the structures (e-g) of Fig. [2] DFT and MCQDPT results 
are represented as black and colored histograms, respectively, according to the labeling system of Fig [2] The color code (for 
MCQDPT results only) is red for majority and blue for minority sites; occasionally, both the singlet (full color) and the triplet 
(shaded) spin manifolds have been considered, as indicated. See text for details. 



4.00 
3.50 
3.00 



^ = 2 



^=1 



2.5) 
2.0) 



1.00 
0.50 



O.OC 



:■! 



E sites 



G sites 



00 1.00 2.00 3.00 4.00 5. 
Np. 



4.00 
3.50 
3.00 
2.50 


CVJ 

< 1 < 1 , 1 , 


1 ' 1 ' 

E sites 1 


2.00 


t 1/ 


£-0 _ 


f 








^ 


■ 


■ 


1.00 
0.50 






_ 


■ ■ 


G sit 


es 


0.0g_ 




1 


, , , , - 


00 


0.05 


0.10 


0.15 O.J 



00 

Np. 



Figure 8. Binding energies as functions of the site normalized 
populations of the substrate HOMO. Black, red and green 
symbols for structures (e-g), respectively. Also indicated the 
hyper coordination number of the edge sites. Left and right 
panels for majority and minority sites, respectively. 
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adsorption site. This is a majority site with an effective 
coordination number Z = 1, which would correspond to 
an additional type of site, "Z)". Furthermore, five E sites 
exist in H-coronene with ^ = 1 having a large binding en- 
ergy. Analogous results holds for pyrene, see Rasmussen 
et d^^. 



Figure 9. Hydrogen adsorption paths on pyrene (a-b) and 
coronene (c-d), on the left for an edge site and on the right 
for a graphitic site. Black and green symbols for DFT and 
MCQDPT results. Lines are spline interpolation to guide the 
eyes. 



E. Adsorption profiles 

We now look at the full energy profiles (minimum en- 
ergy paths) for a H atom adsorption, focusing on few 
illustrative cases. We show in particular that the argu- 
ments used so far for the adsorption energies equally ap- 
ply to the energy harriers for the H atom sticking. Thus, 
the energy ordering rules drawn in the previous sections 
not only determines the thermodynamics but also the ki- 
netics of the hydrogenation process. 



Hydrogen atom binding is an activated process with an 
energy barrier which typically prevents adsorption un- 
der room temperature conditions'^ ^^. For instance, in 
graphite (graphene) the barrier is '^0.2 eV high and this 
prevented for some time observation of a chemisorbed 
hydrogen phase. This barrier is typically linearly related 
to the binding energy itself'^, in accordance with the 
general finding (known as Br0nsted-Evans-Polayni rule) 
that a larger reaction exothermicity is accompanied by a 
lower energy barrier. The same applies here, as is shown 
for the cases of pyrene and coronene reported in FigJ9| 
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Figure 10. Adsorption curves for a second H atom in the 
(graphitic) ortho- (a), meta- (b) and para- (c) positions with 
respect to a first H atom, as indicated in the insets. Grey and 
colored symbols for DFT and MCQDPT results, respectively, 
and lines are spline interpolation to the data for guiding the 
eyes. Panel (d) shows the spin-density of the H-coronene sub- 
strate with a H atom adsorbed on a C site. 



for both an E and a G site. Such curves have been ob- 
tained by fixing the CH distance at the desired value and 
performing a full structural relaxation of the remaining 
degrees of freedom at the DFT-B3LYP level of theory. 
As is evident from the figure, a larger binding energy 
reflects a smaller adsorption barrier, which can be even 
almost vanishing when H binding occurs at an edge site. 
Similar results hold for all the paths considered in this 
work, i.e. for H atom adsorption on most of the sites 
considered in FigJ2] As already noticed above this find- 
ing suggests that the edges of realistic samples could be 
active sites where hydrogenation starts and propagates 
into the bulk: addition of H atoms to E sites modifies 
the sublattice imbalance and at the same time effectively 
converts a number of F (G) sites into E (F) sites. 



function calculations to give a larger binding energy than 
DFT for the graphitic sites, see e.g. the right panel of 
Fig. |4] This is particularly evident for the G site of 
the coronene molecule: Fig. |9] (d) shows that binding 
to this site is '^0.2 eV stronger when computed at the 
MCQDPT than at the DFT level of theory, and that a 
corresponding trend is found for the barrier. However, 
given the limited number of active electrons that could 
be consistently included in the wavefunction calculations 
we doubt that this discrepancy is a manifestation of a 
true physical effect. This is made more evident in Fig. 
[Toj where the adsorption paths for a second H atom onto 
the ortho- ^ meta- and para- position to the first G sites 
are displayed for the two different levels of theory. We 
chose to focus on this system because of the role it played 
as a cluster model for graphene (graphite) since Jeloaca 
and Sidis^^ used it to investigate H atom adsorption on 
the graphitic sites. As is clear from Fig. [To] the above 
discrepancy doubles when adsorption proceeds in para- 
but vanishes for the ortho- site, thereby suggesting that 
the "extension" of the structure may be a source of error 
in the MCQDPT calculations. Notice that the wavefunc- 
tion calculations are always two-state MCQDPT calcu- 
lations, in order to correctly handle the barrier region, 
and included in some cases a level shift correction to get 
rid of the intruder state problem. 

Finally, we performed few additional calculations of 
the binding energies with a very different implementa- 
tion of the DFT-GGA theory, namely a F-point, periodic 
plane- wave calculation using a pure GGA functional as 
described in Section [Till We find for coronene 1.42 and 
0.67 eV for adsorption on the E and the G site, respec- 
tively, which compare very well with the values obtained 
with the hybrid B3LYP functional, namely 1.42 and 0.61 
eV. The same holds for the adsorption of a second atom 



on the same sites considered in Fig. 10 we obtain 2.04, 
0.70 and 1.82 eV for the ortho- ^ meta- and para- graphitic 
sites, to be compared with 2.05, 0.69 and 1.65 eV. Notice 
that also in this case the larger discrepancy occurs at the 
para- position, which might signal the need of additional 
care in the correlation problem. 



V. SUMMARY AND CONCLUSIONS 



F. Correlation level 

Finally, we focus on some technical aspects concern- 
ing the treatment of electron correlation. Though not 
emphasized so far, the results of the DFT-B3LYP cal- 
culations have been shown in parallel to the results of 
more accurate, though more expensive, MCQDPT cal- 
culations (see Section III) which we performed on the 
DFT-optimized structures. As is evident from Fig.s 4|7 
and |9] the two sets of data agree well with each other, 
the discrepancies being at most few tenths of eV in few 
cases. No general trend is found in the comparison, ex- 
cept maybe for a general tendency of the correlated wave- 



We considered atomic hydrogen adsorption on a num- 
ber of small graphenic structures (PAH molecules) in or- 
der to investigate the enhanced reactivity of the edge sites 
already observed by several authors. To this end, we se- 
lected only small structures to prevent the formation of 
radical species at the edge, as it occurs with the forma- 
tion of zero-energy states at the edges of large zig-zag 
nanoribbons. Surprisingly, we found that some edge lo- 
calization always occurs as a consequence of the reduced 
coordination of E sites which translates into a lower on- 
site energy in a renormalized lattice. Further localization 
occurs when E sites are highly coordinated in the renor- 
malized lattice, as measured by a "hypercoordination" 



number ^. We found a very good correlation between the 
binding (barrier) energies and the coordination and hy- 
percoordination numbers: the most favoured sites for H 
atom adsorption (but hkely for adsorption of any mono- 
valent species used to form covalent bonds with carbon) 
are those showing the lowest coordination and the largest 
hypercoordination numbers {E' sites in Fig. [2]). We also 
found, similarly to graphene, that further enhancement 
of the reactivity of specific lattice positions may arise 
from the same topological frustration which gives rise to 
midgap states, i.e. that occurring when the maximal 
set of non-adjacent sites exceeds half the total number 
of sites. In this case a preference towards the maximal 
set of non- adjacent sites adds to the above preference for 
coordination and high hypercoordination. 

We obtained these results in small (subnanometer- 
sized) graphene structures, but they are expected to hold 
for more complex structures. For instance, hydrogena- 
tion is known to occurs much more easily on a zig-zag 
than on an armchair edge of large area graphene: May 
et a/.E3, for instance, extrapolated DFT values computed 
on finite size graphenes towards the infinite size limit and 
obtained 2.86±0.15 eV for the zig-zag edge and 1.74±0.11 
eV for the armchair one. This is consistent with the 
"rules" found here. Indeed, both edges have F and E 
sites, but only zig-zag E sites can be fully hypercoordi- 
nated: ^ = 2 in this case, to be compared with ^ = for 
the E sites of an armchair edge. 

Beside their simplicity, one of the main advantage of 
the derived rules is that they are based on local consid- 
erations which hold irrespective of the global electronic 
properties of the carbon nanostructure under study. As 
a consequence, our findings suggest that exposing ar- 
bitrarly shaped graphene dots to controlled amount of 
atomic hydrogen {e.g. under cold plasma conditions) hy- 
drogenation starts from the edges and propagates into the 
bulk in a much more efficient way than expected solely 
on the basis of the bulk adsorption energetics. 
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